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Abstract. We consider a non-equilibrium three-state model whose dynamics is 
Markovian and displays the same symmetry as the three-state Potts model, i.e., the 
transition rates are invariant under the permutation of the states. Unlike the Potts 
model, detailed balance is in general not satisfied. The aging and the stationary 
properties of the model defined on a square lattice are obtained by means of large- 
scale Monte Carlo simulations. We show that the phase diagram presents a critical 
line, belonging to the three-state Potts universality class, that ends at a point whose 
universality class is that of the voter model. Aging is considered on the critical line, 
at the voter point and in the ferromagnetic phase. 
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1. Introduction 

The critical behavior of equilibrium systems can be classified according to universality 
classes that are characterized by a small number of features such as symmetries, 
conservation laws and dimensionality [ij. The universality classes of non-equilibrium 
systems without detailed balance are also characterized by these features and possibly 
others [21 [3l HJ |5]. It seems that any universality class that encompasses systems in 
equilibrium includes also non-equilibrium systems although the reverse is not true. 
The universality class of directed percolation for instance does not include systems 
in equilibrium. In other words, some classes of universality comprehend equilibrium as 
well as non-equilibrium systems and in this case we may say that the lack of detailed 
balance is irrelevant. Other classes are constituted only by non-equilibrium systems and 
the lack of detailed balance is a relevant property. 

Here we study a two-parameter non-equilibrium model that belongs to the 
universality class of the three-state Potts model. As happens to the latter, the non- 
equilibrium model studied here displays a continuous phase transition from an ordered 
to a disordered phase. As one varies the parameters along the critical line separating 
the ordered and disordered phases, the model leaves the Potts universality class and 
displays a crossover to the non-equilibrium universality class of the voter model at the 
terminus of the critical line. The model is studied on a regular square lattice and the 
transition probabilities depend only on the four neighbors of a given site. We set up the 
most general model obeying the same symmetry as the equilibrium Potts model. In this 
case the symmetry of the model is identified as the symmetry of the transition rates in 
contrast to equilibrium models for which the symmetry of the model is identified as the 
symmetry of the Hamiltonian. The most generic model has five parameters but here we 
focus on a certain two-dimensional section of the space of parameters. 

The model is presented in section [21 The critical line and the critical exponents 
in the stationary state are determined by Monte Carlo simulations and the results are 
presented in the section [3l Finally, the aging of the model when quenched from the 
high-temperature phase is studied in section HI In order to define a response function, 
the coupling to an external field is implemented by an appropriate modification of the 
transition rate as detailed in the Appendix. 

2. Model 

We consider a general non-equilibrium three-state model on a square lattice with full 
symmetry Z^. Any permutation of the three states leaves the transition probabilities 
invariant. On a square lattice, in which each site has four nearest neighbor sites, the 
most general model has five parameters as shown in table [2l Several non-equilibrium 
models as well as the stochastic equilibrium Potts model [6] are particular cases of this 
general model. 

Among them we found the non-equilibrium majority Potts model [71 [8] defined in 
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such a way that a site takes the state of the majority of the neighboring sites with 
probabihty p and any one of the other state with probabihty (1 — p)/2. With our 
definition, it is recovered when 

Pi = P2 = P3 = P4 = 1 - 2^5 = p. (1) 

This model does not obey detailed balance and therefore displays irreversibility at 
the stationary state. It has been shown [7j that it presents a phase transition from 
a disordered (paramagnetic) to an ordered state (ferromagnetic), for sufficiently large 
values of p, which is similar to the phase transition in the ordinary equilibrium three- 
state Potts model. Both models are in the same universality class with the critical 
exponents u = 5/6 and (3 = 2/15. 

The equilibrium Potts model [6] with a Glauber dynamics is recovered when the 
parameters are such that detailed balance is fulfilled. This occurs when 

p. = (2) 



r.3 

l + r + r-^' ' " 1 + r + r'-^ 



2r2 r 
1 + 2r^ 2 + r 

where r = e^'^ and J is the nearest-neighbor exchange coupling of the Potts model. 
Another particular case is the so called linear voter model obtained when we set 

Pi = ^(l + 27), (5) 
P2 = ^(l + ^7), p5 = ^(l-^7), (6) 

P3 = ^(2 + 7), P4 = ^(l + ^7). (7) 

This model is particularly interesting because it can be solved exactly [16]. It displays 
a disordered phase except at 7 = 1 corresponding in our notations to 

Pi = l, P2 = 3/4, P5 = l/4, P3 = l, Pa = 1/2. (8) 

In this particular case it is called the voter Potts model and undergoes a first-order phase 
transition with a magnetization jump but divergences of susceptibility and correlation 
length associated to critical exponents 7 = 1 and u = 1/2. 

The general model includes also models with absorbing states. Due to the Potts 
symmetry these models display three absorbing states represented by a lattice full of 
each of the three states. The class of absorbing models occurs when pi = 1. A particular 
case of this class is the Potts voter model already mentioned and defined by Eq. ([H]). 

Here we are interested in a subclass of models that interpolates between the majority 
vote model, defined by ([I]), and the voter model (IHl). This subclass of models is defined 
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Table 1. Transition probabilities of the general non-equilibrium three-state model on 
a square lattice. The first column shows the states of the four nearest neighbor sites. 
The first row shows the final states. The parameters pi, p2, Ps, Pi, and p^ are restricted 
to the interval [0,1]. The other transition probabilities are obtained by permutation of 
the states. 





1 


2 


3 


1 1 1 1 


Pi 


(l-Pi)/2 


(1 


-Pi)/^ 


1112 


P2 


P5 


1 - 


-P2-P5 


1122 




P3/2 


1 - 


-P3 


1123 


P4 


(l-p4)/2 


(1 


-Pa)/2 



by a set of two parameters a and b that corresponds to a section of the full space spanned 
by the parameters {pi}. The relation between a, b and {pi} is given by 

Pi = b, p3 = b, p4 = a, (9) 

P2 = bi^l + ^y p,= ^-{l-a) + ^-{l-by. (10) 

The phase diagram, to be presented shortly, is therefore a bidimensional section of the 
full five-dimensional space. The majority Potts model is recovered when a = b. When 
6 = 1, we get a model with absorbing state with the particular case of the voter Potts 
model occurring at 6 = 1 and a = 1/2. 

The construction that is presented above is easily generalized to any number of 
states q of the Potts model. Such a construction has already been proposed in the case 
g = 2, i.e., with the symmetry of the Ising model |9l|10]. Among the models appearing as 
special case in the q = 2 construction, voter and linear models were studied analytically 
[TT] and the Glauber-Ising model was studied by Monte Carlo simulations [9| [TOl [T2| [13] . 
Aging was also considered [H]. 

3. Critical properties in the stationary state 

In this section, we give numerical evidences that the phase diagram in the a — b plane 
presents a critical line separating a ferromagnetic phase from a paramagnetic one and 
ending at the point a = 1/2 and 6=1. Excluding this special point, the universality 
class of all points along the critical line is that of the three-state Potts model. Along the 
line b = 1, the system undergoes a dynamical absorbant-active phase transition. The 
intersection of the two lines, at a = 1/2 and 6=1, belongs to the voter universality 
class. 

3.1. Numerical methodology 

The phase diagram of the model was investigated by means of Monte Carlo simulations 
on a square lattices with N = Li^ sites and periodic boundary conditions. The system 
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is initially prepared in the ferromagnetic state, rithemi. Monte Carlo steps (MCS) are 
performed to let the system reach a stationary state. Then niter. MCS are made to 
sample the stationary probability distribution. The order parameter, abusively called 
magnetization in the following, is measured for a given spin configuration a = {cTj} as 

m{a) = ^— ^ , (11) 

where g = 3 is the number of Potts states and Pmax is the density of spin in the majority 
state. If we denote by (Xmax the majority state (among the q possible states) then 

Pm£«((7) ^J^Yl <^^i,^max- (12) 

i 

Of particular interest for the determination of tiie critical line are the susceptibility 

X = N{{m')-{mf), (13) 
and the Binder cumulant U 

C'=l-#4- (14) 

In addition to the susceptibility and Binder cumulant, we also evaluated the 
correlation length. For a finite system, the correlation length ^ is finite so that the 
spatial correlation function decays exponentially, 

C(f, f) = {m{r)m{f^)) - (m)^ + ^e-'^^-^'l/^ (15) 

where m(rj) = (q'^o-i.o-max ~ 1) is the local magnetization on site i. The CPU time 

required to compute the full spatial correlation function being prohibitive, we restricted 
ourselves to the first coefficient of its Fourier transform 

C(^) = ^ VC(0,f)e^^-^"'~ — - — . (16) 

which is efficiently computed as 

C{k) ^ {m{k)m{-k)) ^ {\m{k)f), (17) 

— * 

where m{k) is the Fourier transform of the magnetization 
m(k) — — ^y^m(r)e 

Finally, the correlation length is estimated as 

(19) 

where kmm = (27r/L,0) is the smallest wavevector on the lattice in order to minimize 
the terms ak"^ + . . . and to be able to neglect them. In the paramagnetic phase, 
C(0) = (m^) while in the ferromagnetic phase, one has to remove the unconnected 
term, i.e. C(0) = (m^) — (m)^. Practically, we will use always the second definition 
unless it leads to a negative value which often occurs in the paramagnetic phase because 



1 — ik.f (^;Lg) 
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Figure 1. Auto-correlation time r as a function of b. The different peaks corresponds 
to distinct values of a. The different colors correspond to lattice sizes L — 24, 32, 48, 
64 and for five values of a only L = 90 and 128. 



of the numerical instability of the difference (m^) — (m)^ between two small and noisy 
values. 

In order to adjust the number of iterations ntherm. and nuer. and to calculate error 
bars, we measured autocorrelation functions of defined as: 

C{t, s) = m\t)m\s) - {mY. (20) 

Assuming that these autocorrelation functions decay exponentially, C{t,s) ~ e^'*^'*'/'^, 
the so-called integrated relaxation time is estimated as 

]^ 'lite 



T 



'"Iter p 



(21) 



'^iter t,s=l 

ntherm. and riiter. are always chosen larger than 20r and lO^'r respectively. Since r gives 
a measure of the number of MCS required to obtain a spin configuration uncorrelated 
from the initial one, the number of effective uncorrelated MCS is of order N/r. As a 
consequence, the error bars can be computed for an observable X as 



AX 



N 



(22) 



The error bars only measure the statistical fluctuations of the data. To detect and 
potentially avoid systematic deviations due to metastable or absorbing states, the 
simulation is repeated n^onf = 20 times for all lattice sizes and the observables are 
averaged over these configurations. In the following, we will denote for simplicity {X) 
this double average. 

The difficulty of these simulations is that the relaxation time r is expected to grow 
very fast with the lattice size, i.e. as r ~ where the dynamical exponent is of order 
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2; ~ 2 (the dynamical exponent will be determined by finite-size scaling at the end of 
the next section) since the dynamics is always local (see figured]). At the largest lattice 
size L = 128, up to nconf = 20 simulations of niter. = 7.3 x 10'' MCS each have been 
performed. Despite the computational effort devoted to these simulations, we cannot 
hope to reach the same accuracy as equilibrium simulations for the three-state Potts 
model for which a cluster algorithm is available (which is not the case for our more 
general model for which detailed balance does not hold). For our largest lattice size, 
L = 128, the auto-correlation time reaches the maximum value of t 7231 on the 
critical line belonging to the q = 3-Potts universality class while it is only r ~ 25 for 
the Potts model with the Swendsen-Wang algorithm |15] . 

3.2. Determination of the critical line 

The critical line is determined from the maximum of the susceptibility. We have 
considered twenty values of the parameter a between a = 1/2 and a = 1 and 
lattice sizes L = 16, 20, 24, 32, 48 and 64. For five specific values of a, we made 
additional simulations for lattice sizes L = 90 and 128. For each value of a and L, we 
have performed Monte Carlo simulations for several values of b as shown in figure |2j 
Additional simulations were performed around the maxima of x- Iii the neighborhood 
of these maxima, we interpolated x{b) with a polynomial of degree two to improve the 
accuracy on the location of the maximum. For each lattice size L and each value a, we 
thus defined a pseudo-critical parameter bc{a, L) as shown in figure El These values are 
finally extrapolated in the thermodynamic limit L — )■ +00 with a polynomial of degree 
two in 1/L. The final critical line bc{a) is presented in figure El 

The same analysis has been applied to the correlation length ^ and to the correlation 
time T. However, since the analysis based on the susceptibility x turned out to be the 
most precise, we did not perform further simulations to refine the location of the maxima 
of ^ and r. 

We have also determine the critical line from the crossing of the Binder cumulant 
U at two different lattice sizes (figure H]). For each value of a and for each pair of 
successive lattice sizes, we used a linear interpolation to get a trial estimation of the 
crossing position. Next, we used five points around the preliminary estimation to 
interpolate a polynomial of degree two from which we get a more precise estimation. 
Note that for a = 1/2, the curves do not cross because be = I. Despite the fact 
that the crossing point does not depend significantly on the lattice size, a more precise 
determination of be is forbidden by statistical fluctuations, which are greater than those 
of the susceptibility and thus give a less accurate extrapolation. In the following, we 
will use the determination of the critical line given by the susceptibility. 

3.3. Critical behavior along the transition line 

We have used finite-size scaling analysis to determine the critical exponents. Figure [51 
shows the maximum of the susceptibility x as a function of the lattice size L for several 
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Figure 2. Susceptibility x as a function of b. The different peaks corresponds to 
distinct values of a. The different colors correspond to lattice sizes L = 24, 32, 48, 64 
and for five values of a only L = 90 and 128. 




Figure 3. Critical line determined from the maximum of the susceptibility x- Different 
lines correspond to distinct lattice sizes and to the thermodynamic limit extrapolation. 
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Figure 4. Critical line determined from the erossing of the Binder cumulant U for 
distinct values of the lattice size. 

values of the parameter a. According to finite size scaling, the maximum of x behaves 
as 

X ~ L-r/r (23) 

We see from figure E] that the case a = 1/2 differs significantly from the other cases. 
The ratio 7/z/ is estimated from the slope of the log- log plot of x versus L. The values 
obtained by this procedure are plotted as a function of a as shown in figure El For 
a = 1/2 we get an exponent compatible with the value 7/z/ = 2 of the voter model [l6j. 
For a > 1/2, the effective critical exponent displays a plateau around 1.78(3) along the 
critical hne, i.e., still slightly above the value 7/z/ ~ 1.733 expected for the three-state 
Potts model. The effective exponent decreases when removing the smallest lattices but 
with the lattice sizes that we have been able to considered, the exact value 7/z/ = 1.733 
is not reached yet. 

At the point a = 1/2, corresponding to the voter model, the susceptibility presents 
a logarithm correction to the dominant behavior [T6] . 

X~[(l-&)ln(l-6)]-\ (24) 

In figure [7] we have plotted x versus [(1 — b) ln(l — b)]~^ for several values of the lattice 
size L. The deviation from a linear behavior is due to finite size effects. Indeed, the 
curves become linear over a wider range of values of [(1 — b) ln(l — b)]~^ as one increases 
L. 

To determine the exponent z/, we have used the position bc{a, L) of the maximum 
of the susceptibility and the value bc{a, 00) obtained by extrapolation. The finite-size 
deviation is expected to behave as 

b^{a,L) -bc{a, 00) ^ L-^^r (25) 
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Figure 5. Maximum of tlie susceptibility x as a function of tlic lattice size L. Different 
curves correspond to distinct values of the parameter a. 
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Figure 6. Effective exponent ^/v obtained by interpolation of the susceptibility as a 
function of a. The dotted lines correspond to ^/v = 2 (voter model) and to ^ /v ~ 1.733 
(three-state Potts model). The different symbols correspond to exponents obtained by 
power-law interpolation when taking into account lattice sizes L > Lmin only. 
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Figure 8. Exponent v obtained by interpolation as a function of a. The dotted line 
correspond to = 5/6 (three-state Potts model). The different symbols correspond 
to exponents obtained by power-law interpolation when taking into account only the 
lattice sizes L > Lmin- 
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Figure 9. Scaling of tlie Binder cumulant U in the critical line. All data points for 
lattices sizes L — 24, 32, 48, 64, 90 and 128 and for values of a above 0.7 were used. The 
data for a ~ 0.57 is also plotted (blue triangles) to show the cross-over with the voter 
point. The two branches, above and below the dashed horizontal line, correspond to 
high and low-temperature phases. 



The exponent u can then be obtained from the log-log plot of the deviation versus L. 
Our numerical estimates are shown in figure El The average value 1/z/ = 1.10(13) leads 
to u = 0.91(10) slightly larger, but not incompatible with the value u = 5/6 of the 
three-state Potts model. 

The scaling of the Binder cumulant is presented on figure M For parameters a 
smaller than a ~ 0.7, a strong cross-over is observed: the data do not fall on a single 
curve for different values of a. For a > 0.7, the cross-over effects become smaller than 
the error bars and all data fall on a unique curve, providing further evidence that the 
universality class is the same along the critical line, at least for a > 0.7. From the curve, 
we can see that the universal value of the Binder cumulant isU = 0.615(5) in agreement 
with [8]. 

On figure [TOl the scaling function J-'^ of the correlation length ^ such that 

e~LJ-5((6-6,(a))L^), (26) 

is plotted for values of a larger than 0.7. Again a strong cross-over is observed below 
a ~ 0.7. In the finite-size regime {b — bc{a))L'^ < 0.1, the scaling function displays a 
plateau at a value ^/L that is expected to be universal. The values of a considered in 
the plot leads to a compatible value of ^/L. Note that our value of C,/L differs from that 
given in ref. [15] for the 3-state Potts model because of a different correlation function 
used for the definition of the correlation length ^. 
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Figure 10. Scaling of ttie correlation length ^ on the critical line. All data points for 
lattices sizes L = 24,32,48,64,90 and 128 and for twelve values of a above 0.7 were 
used. The points on the second branch at the right of the figure correspond to data 
in the paramagnetic phase for which for numerical reasons (see previous section), the 
definition of ^ has to be chosen differently. 



To determine the dynamical exponent z we consider the maximum of the auto- 
correlation time r. According to finite size scaling this quantity behaves as 



Figure [TT] shows the log-log plot of this quantity versus L for several values of a. The 
dynamical exponent z, estimated from the slope of these curves, is shown in figure [T2j 
We obtain an exponent z ^ 2.25(6) compatible with the recent values, around z ~ 2.197, 
obtained for the three-state Potts model pTl [T8] . The value z = 2 of the voter model 
is reproduced by the data at a = 1/2 for large lattice sizes only. A strong cross-over is 
observed for values of a close to 1/2. The effective exponents z smaller than 2 close to 
a = 1/2 are artificial: these too small values are due to the fact that the auto-correlation 
times display first a regime in the voter universality class and then a second one in the 
Potts universality class. But for the same value of L the auto-correlation time is larger 
at the voter point than in the critical line so that in the cross-over region, r has to first 
reduce its growth to join the second regime (see figure [TTl) . 

4. Aging 

4.I. Numerical methodology 



(27) 



The system is now considered out of the stationary regime, that is, in the transient 
regime occurring before it reaches the stationary state. It is initially prepared in a 
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Figure 12. Dynamical exponent z obtained by interpolation of the auto-correlation 
time T. The different symbols correspond to exponents obtained by power-law 
interpolation with taking into account lattice size L > Lmin only. 
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random state, i.e., in the paramagnetic phase. At time t = 0, it is "quenched" to 
some values of a and b. A number t/ of Monte Carlo steps are then performed. The 
numerical experiment is repeated riconi times. In this section, X{t) means an average of 
the observable X at time t over all these nconf histories. 

Motivated by the fact that it is possible to couple a magnetic field to any of the 
q spin states and thus to define q possible response functions, we will now define q 
spin-spin autocorrelation functions as 



C^{i;t,s) = m^{i;t)m,^{i;s), (28) 

where m„{i;t) = [q6^-(^t).a ~ ~ 1) is the local magnetization at time t on site i. 

Since translation invariance is expected, we can average over the lattice to get a more 
stable estimate: 

CAt,s) = j^Y.CAr,t,s). (29) 

i 

Finally, all Potts states being equivalent, one can define a spin-spin autocorrelation 
function as 

C{t,s) = j2CAt,s). (30) 

cr=l 

One can easily check that these definitions lead simply to 

C{t, ^) = (g ^ 1)2 E ^Mt)Ms) - • (31) 

For uncorrelated spins at times t and s, this autocorrelation function vanishes. 

A magnetic field h is coupled to a particular value ah among the q Potts states by 
modifying the transition rates by a multiplicative factor [161 IlHl ED] : 

w{a, al) = wo{ai ^ a',)e-''^-^-'^^ . (32) 

The magnetic field is branched only during one MCS, at time s. The linear response of 
magnetization to this field is defined as 



where again m„{i]t) = (q6„.(^t),(T ~ 1)/(q' ^ 1) is the local magnetization. The numerical 
method used to calculate this response is presented in the Appendix. As for the 
correlation function, the response can be averaged over the lattice and then the different 
responses can be added since all states are equivalent: 

Rit,^) = ^Il T.R^Sr,t,s). (34) 

Finally, the Fluctuation-Dissipation Ratio (FDR) is defined as 

Xit s) - ^ (35) 

^ ' ' d,C{t, s) - C{t, s + 1)- C{t, s) ■ ^ ^ 
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We have also considered global quantities: the magnetization-magnetization 
correlation function 

CMit, s) = M{t)M{s) = ^ E rn.(^;t)m^(j;s), (36) 
and the response to global magnetic field 

as well as the associated Fluctuation-Dissipation Ratio. 



4-.2. Aging at the voter point 

We consider first a quench at the point a = 1/2 and 6 = 1 belonging to the voter 
universality class. As will be shown in the following, the theoretical predictions given 
in Ref. [16] are well reproduced by the numerical data. The spin-spin two-time auto- 
correlation function C(t, s) is expected to behave as 

We used a lattice size L = 192, a final time tf = 1000 and observables have been averaged 
over riconf = 10^ histories of the system. As shown on figure [13], the data indeed tend 
toward the theoretical prediction when t and s are sufficiently large. The exponent 
X/z is thus 1. The Fluctuation-Dissipation Ratio X{t,s) is presented on figure Its 
asymptotic limit is compatible with the theoretical prediction X^o = 1/2. Interestingly, 
the FDR associated to the global magnetization reaches a plateau at the value 
almost immediately. Even though, this quantity is noisier, it can be used to estimate 
at very short times at the price of a higher statistics. 



4.3. Aging on the critical line 

We studied the aging at several points on the critical line. Because of the cross-over 
with the voter fixed point, several lattices sizes, L = 64, L = 128 and L = 256, were 
considered. 

The spin-spin correlation function is expected to behave asymptotically as [211 122] 

/t \ ^^/^ 

C{t,s)^s-'^^(^-j , (39) 

where Oc = 2(5 /vz. The scaling form is tested by plotting the scaling function s°'''C{t, s) 
with respect to t/s. The collapse of the curves for different waiting times is obtained 
(see figure [T5D but the exponent Oc leading to this collapse tends to grow along the 
critical line: from the expected value 2(5 /uz ~ 0.11 close to the voter point to 0.2 
at a = 1. This indicates that for the largest values of s, the asymptotic regime 
C{s,s) ~ s"""" is not reached yet and corrections to this dominant behavior cannot 
be neglected. Nevertheless, the resulting scaling function s"=C(t, s) appears very similar 
along the critical line and an algebraic decay with t is observed for all waiting times. 
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Figure 13. Scaling function lnsC(i,s) versus t/s— 1 at tlie point a — 1/2 and h—\ 
belonging to the voter universality class. The dashed line is '^/{t/s — 1). 




Figure 14. Fluctuation-Dissipation Ratio X{t, s) versus s/t at the point a = 1/2 and 
h = 1 belonging to the voter universality class. In the inset, the FDR associated to the 
global magnetization is plotted. 
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Figure 15. Scaling function s°'''C{t, s) with respect to t/s for two points on the critical 
line (a — 0.63158 on the left and a = 1 on the right). The different curves correspond 
to different waiting times s. The lattice size is L = 256. 
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Figure 16. Effective exponent — A/z as a function of the smallest bound t/s of the 
interpolation window of C{t,s). From left to right, a = 0.63158, 0.89474 and a = 1. 
The different curves correspond to different waiting times s. The dashed line is the 
power-law regression of the data. 



This provides further evidence of a unique universahty class along the critical line. The 
exponent X/z is obtained by power-law interpolation of C{t,s) with the time t. For 
each value of the waiting time s, the data are interpolated within a sliding window 
[t/s,tf/s]. The effective exponents are represented as a function of s/t in figure fT6l 
The asymptotic value is read at the intercept with the ?/-axis while the small slope 
that can be observed comes from the corrections to the behavior {t/s)~^^^. We finally 
extrapolated the effective exponents and obtained X/z ~ 0.806(5), 0.817(4), 0.820(4) 
and 0.818(3) for a = 0.63158, 0.76316, 0.89474, and 1, respectively. These values should 
be compared with the values 0.828(2) [17] and 0.844(19) [23] obtained in the case of the 
pure three-state Potts model. 

The same exponents can be obtained from the scaling of the response function 
R{t,s). In contradistinction to the correlation function, the collapse of the scaling 
function s^^""'R{t,s) is obtained with the expected value Oc = 2(5 /vz all along the 
critical line (see figure [T7|l . The response function seems to be affected by weaker 
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Figure 17. Scaling function s^^"-''R{t,s) witli respect to t/s for two points on the 
critical line (a — 0.63158 on the left and a = 1 on the right). The different curves 
correspond to different waiting times s. The lattice size is L — 256. 




Figure 18. Fluctuation-Dissipation Ratio X{t,s) versus s/t for a — 0.63158 (left) 
and a — 1 (right). The different curves correspond to different waiting times s. 



corrections than the correlation function. However, the response is noisier than C(t, s) 
and the effective exponents \/z obtained from a power-law interpolation of the decay 
of R{t, s) with t are compatible but much less accurate than the exponents computed 
from C(t, s). 

From the response function and the derivative of the correlation function, the 
Fluctuation-Dissipation Ratio X(t, s) was computed. The asymptotic value is 
expected to be universal. As the parameter a is increased, larger and larger fluctuations 
are observed (figure [T8|) . The asymptotic value seems to be a bit lower, at least for the 
small waiting times, than the value X^o — 0.406 for the three-state Potts model [23]. It 
remains nevertheless always compatible within statistical fluctuations. 

4.4- ^ging in the ferromagnetic phase 

The system is now prepared in a random state and quenched in the ferromagnetic phase 
at the point of parameters a = 0.9 and h = 0.95. The auto-correlation exponent X/z 
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Figure 19. Effective exponent X/z obtained by a power-law interpolation of the 
autocorrelation function in a shrinking window [t;tf] and plotted with respect to 
(s/ty^^. The dashed line is the linear interpolation of the s = 320 curve. 



is extracted from the decay of the autocorrelation function C{t,s) ~ M'^{s)t^^^^ . The 
corrections to the dominant behavior appear to be very strong. Like on the critical line, 
we have computed an effective exponent A/2; by power- law interpolation of C{t, s) within 
a shding window [t/s,tf/s]. The asymptotic value is expected to be recovered when t 
goes to tf. It turns out that the effective exponent behaves linearly with (s/t)^/^ for 
sufficiently small values oi s/t (figure [T9|) . This simply means that the autocorrelation 
function behaves as C()f:, s) ~ t^^/^[l + at~^/^ + . . .]. The exponent A/2; is extrapolated by 
a simple linear fit of the data in the linear regime. The value 0.595(15) obtained for the 
three-state Potts model [21] is nicely reproduced for the largest waiting times: we get 
0.605(3) for s = 80, 0.599(3) for s = 160 and 0.601(3) for s = 320. These exponents are 
remarkably stable when the interpolation window is changed (at least when remaining 
in the linear regime). In the ferromagnetic phase, the Fluctuation-Dissipation Ratio 
X(t, s) is expected to decay as to zero. The numerical data confirms this prediction 
(figure [20]) with the predicted value a = 1/2. 



5. Conclusions 



We have shown that the class of non-equilibrium three-state lattice models with Z3 
symmetry that we introduced in this paper displays a phase diagram consisting of a 
critical line belonging in the three-state Potts model universality class and of an ending 
point in the voter class. The simulations for this model are difficult because of the 
absence of cluster algorithm that would reduce the critical slowing down and because 
the critical line cannot be determined by self-duality arguments. We thus obtain less 




Figure 20. Fluctuation-Dissipation Ratio X{t, s) with respect to s/t when the system 
is quenched in the ferromagnetic phase. 



accurate estimates of the critical exponents than for the usual equilibrium three-state 
Potts model. Therefore, we essentially limited ourselves to check the compatibility with 
the theoretical predictions, namely, that non-equilibrium models with the symmetry of 
the Potts model belongs to the universality class of the equilibrium Potts model along 
the critical line. The cross-over with the point in the voter universality class remains 
weak in the Finite-Size Scaling of the averages in the stationary state but turns out 
to be very strong in the dependence of the scaling function with b. In a second step, 
the model was studied in the aging regime after a sudden quench either on the critical 
line or in the ferromagnetic phase. The results are again compatible with either the 
theoretical predictions for the voter model or the numerical estimates for the Potts 
model. Surprisingly, the corrections are stronger in the ferromagnetic phase than on the 
critical line. 
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Appendix: Numerical estimation of the response 

Several different methods for the calculation of the response function without applying 
any magnetic field have been proposed. Our method is based on [25] and was adapted 
to our more general transition rates. For completeness, it is presented in the following. 
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Preliminaries on master equations 

The dynamics is a Markovian process governed by the master equation 

dtVia, t)=J2 [^(^', t)W{a' ^ a) - t)Wia ^ a% (40) 

cr' 

where a denotes the vector whose components are the spin variables {cTj} attached to the 
sites of the lattice. The transition rates W{a ^ a') do not appear in the equation when 
0" = cr' so they are free to take any value. In discrete time (Monte Carlo simulation), 
this master equation is replaced by 

V{(r, t + l)= V{a, t)+Y. [^(^'' ^)^(^' ^ ^) - t)W{(T ^ a')], (41) 

cr' 

which can be written as 

V{a,t + 1) =J2'Pi(^\t)u{a' ^ a), (42) 

cr' 

where we have introduced the quantity 

uj{a' ^a) = [l-Y.Wia^ a")]6^,^, + Wia' ^ a), (43) 

cr" 

that can be interpreted as the conditional probability for the system to undergo a 
transition from a' to a if a 7^ a' or the probability to stay in the same state if a = a'. 
It is easily shown that 

cr 

as required for the probability V{a, t) to be normalized at any time t. 



Response function in the general case 

We are interested in the response to a magnetic field 



where mj(t) is the average magnetization on site i at time t. The response to a global 
magnetic field H = J2ihi acting on all sites is measured as 



dM{t) dmi{t) dmi{t) ^ 

To measure the linear response to a magnetic field hj at time s, the transition rates are 
modified between the interval of time [s; s + 1[. This is obviously done with Glauber 
transition rates by introducing a Zeeman interaction in the Hamiltonian. In the present 
case of non-equilibrium models, that are not defined by a Hamiltonian, the perturbation 
should be introduced into the transition rate. We use here a prescription similar to that 
used for the non-equihbrium models with two states [El EO] and which has been used 
in Ref. [16], namely, 

Wh{a ^ a') = W{a a') g-^. (47) 
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Notice that this choice does not break the detailed balance if it was satisfied by 
W{(7 — 7- cr'). Other choices of perturbation include the one introduced in Ref. [26] . 
The average magnetization on site i and time t is given by 



cr 

= mi{a)r{a, t\a' , s + l)V{a', s + 1) 

o-.cr' 

= rni{(T)V{a,t\a',s + l)io{a" ^ a')V{a",s). 



(48) 



When branching the magnetic field, the transition rate 0Jh{cr' — t- a) has to be replaced 
by 



so that the average magnetization at time t is perturbed by the magnetic field 



TTiiit) = m,{a)V{a, t|a', s + l)uoh{a" a')V{a" , s). 

(T,<T',cr" 

The linear response follows: 

R.,{t,s)= Y m,{a)Via,t\a',s + l)^{a" ^a')Via",s) 

III I in 



(50) 



Y mi(a)P(a,t|a',s + \)\j^m,{a^)W{a' -> a-^b,, ,,nV{a\s 



0-3 



Y m,{a)V{a, t\a', s + l)mj{a')W{(T" (r')V{a", s). 



(51) 



From the definition of u){a" — > a'), the transition rate can be written as 

W{a" a') = Lo{a" ^a')-[l-Y ^ ^s)]'^^',^". (52) 

Introducing this equality into the second line of the expression of the response, one gets 



Rij{t,s)= Y "^i((T)P(a,t|a', s + 1) 

a,a',a" 

+ H ('"T'ji^^s) - mj{a'))W{a' (y3)\K',a"V{a" , s) 

- Y rn,{a)Via,t\a',s + l)mj{a')uj{a" ^ a')V{a",s). 

In the second line, the sum over a" can be performed. It leaves a correlation function 
of observables measured at time t and s + 1. However, in the first line, a transition rate 
u between the spin configurations at time s and s + 1 is missing to do the same thing. 
We thus write 

Rij{t,s)= Y m{(r)V{a,t\a',s + 1) m,j{a') + Yi^j{(^) -'^ji(^')W{(^' ^ 

5.. 



uj(a' — )■ a 



a')V{a\ s) 
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m,{a)V{a, t\a', s + l)mj{(T')u{(T" ^ (t')V{(t", s). 



Since 

the response is finally given by the general expression 



Rij{t, s) = rriiit) /\mf^{s)5^(s),a{s+i) - mj{s + 1) 



where 



^ mj{s) + {mj{(y') - mj{s))W{a{s) a') 
~ W{a{s) a{s)) + 1 - Ea' W{a{s) a') 
In the particular case where the transition rates are chosen such that 

a' 

which is the case in this work, the response reduces to 



(53) 

(54) 
(55) 
(56) 

(57) 

' W{a{s) ^ a{s)) • ^^^^ 

In the case of the Ising model with Glauber transition rates, one recovers the result 
of [25] up to a factor —(5 due to a difference in the definition of the magnetic field. 
Note that is not possible to calculate the response in this way if there exists a spin 
configuration a for which Vr((j — )■ a) = 0. The trick to circumvent this difficulty, is to 
modify the model as 

W\a ^ a') = p6^y + (1 - p)W{a ^ a'), (59) 

The dynamics is then simply slown down by a factor p. 



where 



Rij{t, s) = rriiit) \mf^ {s)5„{^s),a(s+i) - mj{s + 1) 



mAs) + 
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